% tic
if species_lobster1_finfish2==1
   EF_E_F_weight_range=[0 0.6 0.7 0.715 0.73 0.74 0.75 0.76 0.77 0.8 1]; 
elseif species_lobster1_finfish2==2 
%flounder
EF_E_F_weight_range=[0 0.4:0.05:0.5 0.52 0.5225 0.523 0.5275 0.532 0.535:0.005:0.54 1]; %0=al energy
end

Dims=length(find(energy_binary_i==1));
MaxSingleChangeTrials=1e4;
TabuOnOff=1; 

energy_i_EF_E_F=zeros(numpatches,length(EF_E_F_weight_range));
energy_EF_E_F_weight=NaN(1,length(EF_E_F_weight_range));
Profit_equil_wrt_Energy_P1P2_Energy_Fishery_EF_E_F_weight=energy_EF_E_F_weight;
Whale_wrt_Energy_P1P2_Energy_Whale_EF_E_F_weight=energy_EF_E_F_weight;

EF_E_F_weight_index=0;
for EF_E_F_weight=EF_E_F_weight_range;
    EF_E_F_weight_index=EF_E_F_weight_index+1;
    disp(['----------EF search weight = ',num2str(EF_E_F_weight),])
    Fsum=Fsum_opt;
    a=BinaryTabuMonteCarlo(@Energy_Fishery_EFvalue, Dims, MaxSingleChangeTrials,TabuOnOff);
    energy_i_EF_E_F((energy_binary_i==1),EF_E_F_weight_index)=a;
    
    energy_EF_E_F_weight(EF_E_F_weight_index)=sum(energy_profit_i(energy_binary_i==1).*a); 
    disp(['Energy proportional value = ',num2str(energy_EF_E_F_weight./max_energy_profit),])
    
    energy_locs=zeros(numpatches,1);
    energy_locs(energy_binary_i==1)=a;


% % % % %%%%%%%%%%%%%%
% % % % % warning off
% % % % % options_OptFsum=optimset('Display','off');
% % % % % [Fsum_opt maxFleetProfit flag]=fmincon(@OptFsum,x0,[],[],[],[],LB,UB,[],options_OptFsum);
% % % % % warning on
% % % % %%%%%%%%%%%%%%
x0=Fsum_opt; %set guess at previous output
path('\\tsclient\H\Data & Models\npsol;',path); %This enables npsol
prob.guess = [x0];
prob.obj = 'OptFsum';
prob.lo = [LB];
prob.up = [UB];
[maxFleetProfit, Fsum_opt, g, hess, flag, iter]=npsol(prob);  % OPTIMIZATION
%%%%%%%%%%%%
disp(['Fsum_opt = ',num2str(Fsum_opt),'; maxFleetProfit = ',num2str(-maxFleetProfit),'; Exit flag = ',num2str(flag),])
Fsum=Fsum_opt; %use new Fsum in next EF Binary opt
MOPmodel2

Profit_iy=Profit_iy.*energy_binary_i_P1P2_repmat; %ONLY VALUE IN P1 AND P2
Profit_y=sum(Profit_iy,1);
Profit_i=Profit_iy(:,end);
if species_lobster1_finfish2==1 %lobster fishing buffer during construciton phase
    Profit_i_construction=Profit_i.*(energy_locs==0);
    Profit_equil_wrt_Energy_P1P2_Energy_Fishery_EF_E_F_weight(EF_E_F_weight_index)=sum((sum([repmat(Profit_i_construction,1,Construction_time) repmat(Profit_i,1,turbine_lifetime)],1)).*discount_lifetime);
%     Profit_equil_wrt_Energy_P1P2_Energy_Fishery_EF=Profit_equil_wrt_Energy_P1P2_Energy_Fishery_EF/divisor;
elseif species_lobster1_finfish2==2
    Profit_equil_wrt_Energy_P1P2_Energy_Fishery_EF_E_F_weight(EF_E_F_weight_index)=sum(Profit_y(end).*discount_lifetime);
%     Profit_equil_wrt_Energy_P1P2_Energy_Fishery_EF=Profit_equil_wrt_Energy_P1P2_Energy_Fishery_EF/divisor;
end
disp(['Fishery proportional value = ',num2str(Profit_equil_wrt_Energy_P1P2_Energy_Fishery_EF_E_F_weight./Profit_equil_wrt_Energy_P1P2(1,1)),])

Whale_wrt_energy_i=sum(sum([repmat((whalew_density_P1P2.*(energy_locs==0)),1,Construction_time) repmat(whalew_density_P1P2,1,turbine_lifetime)],1).*discount_lifetime); %discounted whale value
Whale_wrt_energy_i=Whale_wrt_energy_i./divisor;
if species_lobster1_finfish2==1
    Whale_wrt_Energy_Fishing=Whale_wrt_energy_i-Fsum*W_F_system*Whale_wrt_energy_i;
elseif species_lobster1_finfish2==2
    Whale_wrt_Energy_Fishing=sum(Whale_wrt_energy_i);
end
Whale_wrt_Energy_P1P2_Energy_Whale_EF_E_F_weight(EF_E_F_weight_index)= Whale_wrt_Energy_Fishing;

end

disp('Done finding energy-fishery EF')
% disp(['----------EF search took ',num2str(toc/60/60),' hours'])


% %%%%%%%%%%%%
% energy_EF_E_F_weight_scaled=energy_EF_E_F_weight./max_energy_profit;
% Profit_wrt_Energy_P1P2_Energy_Fishery_EF_E_F_weight_scaled=Profit_equil_wrt_Energy_P1P2_Energy_Fishery_EF_E_F_weight/Profit_equil_wrt_Energy_P1P2(1,1);
% figure
% plot(Profit_equil_wrt_Energy_P1P2_scaled(:).*100,energy_wrt_Energy_scaled(:).*100,'o')
% set(gcf,'color',[1 1 1]); box off
% ylabel('Energy [%]'); 
% xlabel('Fishery [%]');
% title('P1 and P2')
% hold on
% plot(Profit_wrt_Energy_P1P2_Energy_Fishery_EF_E_F_weight_scaled.*100,energy_EF_E_F_weight_scaled.*100,'ro')
% axis tight